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We have found a way to analyze Edwards’ density of states for static granular packings in the 
special case of round, rigid, frictionless grains assuming constant coordination number. It obtains 
the most entropic density of single grain states, which predicts several observables including the 
distribution of contact forces. We compare these results against empirical data obtained in dy- 
namic simulations of granular packings. The agreement between theory and the empirics is quite 
good, helping validate the use of statistical mechanics methods in granular physics. The differences 
between theory and empirics are mainly due to the variable coordination number, and when the 
empirical data are sorted by that number we obtain several insights that suggest an underlying 
elegance in the density of states. 
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The intriguing behaviors of sand and other granular 
materials are not well understood from a fundamental 
point of view [l] and there is no theory with a pedi- 
gree equivalent to the Navier-Stokes or Maxwell’s equa- 
tions to explain how their state evolves over time, even 
in the most commonly occuring scenarios. Consider- 
ing the ubiquity of granular materials in nature, this is 
quite surprising. Making a new effort to explain their 
physics, Edwards, et al , have proposed that the meth- 
ods of statistical mechanics may be successfully applied 
[2], They hypothized a priori a flat measure in the sta- 
tistical ensemble — that every metastable arrangement of 
grains (a blocked state) is equally probable under certain 
conditions — and that the analysis of this ensemble should 
predict some of the important behaviors. 

There is no theoretical proof of ergodicity for granu- 
lar materials, and so it is an important question whether 
or not the nonlinear, lossy dynamics of granular mate- 
rials might produce significant selection effects, biasing 
the measure so that Edwards’ hypothesis would not be 
correct. Seeking to answer this, a number of empiri- 
cal tests have been performed by computer simulation. 
In these, Edwards’ hypothesis has successfully predicted 
packing behaviors for several idealized models [3] and 
the diffusion-mobility behavior of individual grains when 
a simulated packing is slowly sheared [4]. In both cases 
it appears that the dynamics cause the geometry of the 
model to explore some region of the phase space with 
sufficient ergodicity to justify the flat measure. 

In this Letter we present a different kind of test for Ed- 
wards’ hypothesis. Rather than examining the geometric 
features of the packing, we demonstrate that Edwards’ 
fiat measure correctly predicts the distributions for sin- 
gle grain load states and for contact forces. This pre- 
diction is centrally important to a statistical mechanics 



FIG. 1: (Semilogarithmic) Distribution of granular contact 
force magnitudes Pf (/) (main graph) and Cartesian compo- 
nents P x {fx) (inset). The theory and the empirical discrete 
element model (DEM) data are strikingly in agreement. Be- 
cause the form of Pf(f) reflects the important organizational 
features of quasi-static granular physics, this result implies 
that Edwards’ hypothesis is sufficient to capture that physics. 


theory because a distribution reflects how the ensemble 
is organized and demonstrates whether or not the correct 
physics have been incorporated. In particular, Edwards’ 
hypothesis should predict at least three features in the 
contact force distribution P/(f) as illustrated in Fig. 1: 
the wide tail [5] related to the heterogeneity of stresses in 
a packing ( force chains) [6]; the small peak near the av- 
erage value of force under isotropic conditions [7] related 
to static equilibrium of the grains ( jamming ) [8]; and the 
non-zero probability density at zero force, Pf( 0) > 0, 
related to the tipping of grains (fragility) [9]. The com- 
bination of these three features is unique to the granular 
distribution, not being found in the typical densities of 
states for thermal systems. If Edwards’ hypothesis fails 
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to produce this form then it is doubtful that it could be- 
come the basis for a theory of quasi-static rheology, since 
the tipping or sliding of individual grains depends upon 
the state of their contact forces. 

This Letter will address the question first by analyzing 
the density of states (DOS) from first principles, which 
results in a tranport equation. We will present the distri- 
butions that are the numerical solution to the transport 
equation. We have performed Discrete Element Modeling 
(DEM) [10] for granular packings under conditions sim- 
ilar to those assumed by the theory. The comparison of 
the two suppports the sufficiency of Edwards’ hypothesis 
and provides new insight into the form of the DOS. 

Following Edwards and coworkers [11], we focus on the 
case of amorphous packings of cohesionless, rigid grains 
all having the same coordination number Z that makes 
the packing isostatic [12]. Our case is further idealized 
by using two dimensional round grains with monodis- 
perse grain diameters, omitting gravity and working in 
the thermodynamic limit (infinitely large packings) so 
that the boundary layer may be neglected. We focus on 
the frictionless case so that .Z — 4, and we limit this Let- 
ter to isotropic stress and fabric although the transport 
equation can solve for anisotropic cases, too. The ide- 
alizations may be taken out in future refinements of the 
theory, but this is a good starting point because packings 
of cohesionless, round grains that are perfectly rigid [13] 
and/or frictionless [14, 15] and/or monodisperse [7, 15] 
have been the focus of many empirical studies and are 
known to have force distributions with the same features 
as the less idealized packings. Hence they must be sub- 
ject to the same basic organization in the physics. 

Our avoidance of friction is especially important be- 
cause it is known that 2D round, rigid, frictionless pack- 
ings in dynamic simulations actually do have an aver- 
age coordination (Z) = 4 such that they are isostatic, 
whereas frictional packings are hyperstatic. Isostatic- 
ity means that a blocked state subjected to a particular 
stress condition will correspond to only one microstate of 
contact forces between the grains. Hyperstatic ensembles 
have been recently analyzed [16], in which each blocked 
state corresponds to more than one microstate of con- 
tact forces — a “solution space” with non-zero volume. It 
has been shown that common dynamic modeling tech- 
niques do not sample that solution space uniformly [17], 
so it would be difficult to compare the forces predicted by 
Edwards’ hypothesis with those resulting from dynamic 
modeling in a hyperstatic case. While the issues sur- 
rounding hyperstaticity are clearly important, we avoid 
them here so that we may test Edwards’ hypothesis. 

The goal of the analysis is to combine Edwards’ micro- 
canonical DOS [18] and contact force probability func- 
tional [11] and derive the density of single grain states 


p g (w Xi w y: Qi ,. . . , 84). The Cartesian loads are 

1 V- A 1 . ^ 

w x = -^2fy\cos9 y \, % = -£/ 7 |sin* 7 |. (1) 

7=1 7=1 

where / 7 and Q 1 are the four contact forces and contact 
angles on a grain. In the special case we have selected, 
solving for p g provides everything that can be known 
about the individual grains in the packing. For example, 
it contains the joint distribution of contact forces and 
angles, 
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and the generalized fabric distribution P4e(0i y . . . ,64) 
[19]. 

The analytical method is to count states in Edwards’ 
ensemble and maximize entropy applying the same well- 
known procedure that has been used to derive the Bose 
or Fermi energy distributions [20]. Several innovations 
were needed to apply this procedure to static granular 
packings as discussed in [21]. The result is 

p s (w x ,w y , e u ...,6 4 ) = G(6 U ... A) 
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where 0 is the Heaviside (unit step) function, X x and 
A y are the Lagrange multipliers that scale mechanical 
loading, and G is the Fabric Partition Factor that derives 
from the array of Lagrange multipliers used to conserve 
P40. Eqs. (2) and (3) form a recursion in P fe and p g , the 
“transport” equation, which may be solved numerically 
using p4e and the mechanical loading as inputs. For the 
present this has been solved for the isotropic case using 
the Mean Structure Approximation (MSA), the details 
of which are in Ref. [21]. 

To compare with the theory, we have performed DEM 
simulations of 17,000 two dimensional, round, friction- 
less grains. A portion of a packing is shown in Fig. 2 to 
contrast the spatial disorder of the force network with 
the simplicity of the statistics, discussed below. The 
grain diameters were uniformly distributed from 1.0 to 
1.5 to reduce crystallization. The grains were deposited 
isotropically into a square test cell with hard walls and 
without gravity, and their diameters increased by rescal- 
ing, producing the desired isotropic stress state. The 
grains were allowed to move dynamically until they lo- 
cated and settled into a blocked state. They have a linear 
spring contact law, but staying near the jamming tran- 
sition avoids excessive deformation of the contacts and 
approximates the grain rigidity of the analysis. Data 
from grains in the boundary layer (chosen to be 4 grain 
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FIG. 2: Portion of a Discrete Element Model (DEM) showing 
the disordered packing fabric and propagation of force chains. 
Although disordered, a simple pattern can be discerned in 
the distribution of load anisotropies, s, and total loads, t. 
These functional forms of these distributions, as well as P/(/) 
and P x {f x ), are correctly predicted by the transport equation 
based on Edwards’ hypothesis. 


diameters along each wall) were discarded to reduce the 
boundary effects, which we found will otherwise signifi- 
cantly skew the statistics. 

Fig. 1 shows the DEM results compared to the the- 
ory for P/(f) and for the distribution of Cartesian com- 
ponents of force, P x (fx)- The two inflections in P x (f x ) 
are characteristic of modified Bessel functions of the sec- 
ond kind, and so we see that the theory successfully pre- 
dicts the functional form. In the isotropic case P/(f) can 
be analyticall derived from P x (f x ), and so its functional 
form must be an exponential with a polynomial prefac- 
tor [21, 22]. It appears there is some error in the region 
of weak forces (better seen on a linear plot), mostly at- 
tributable to the theory’s constant Z = 4 for every grain 
as further discussed below. In the DEM there are signif- 
icant populations for Z = 3, 4, and 5. 

To convey the other p g information from this solution 
we might note that w x and w y are not statistically inde- 
pendent and therefore plot their statistics as a joint dis- 
tribution. However, the mean structure theory predicts 
that the change of variables to total load , t = w x +w y , and 
load anisotropy, s = (w x — w y )/t, achieves (approximate) 
statistical independence so that they can be separated 
more meaningfully. The distribution of load anisotropy, 
P s (s), is shown in Fig. 3, demonstrating remarkable 
agreement between DEM and theory. We can fit a func- 
tional form to both, P s (s) = cos(7rs/2) exp(— 8s 2 ). The 
cosine arises in the transport equation due to the struc- 
tural connectedness of the packing, and the Gaussian fac- 
tor due to a grain’s stability space. We segregated the 
DEM’s grains into populations with Z = 3, 4, or 5 and 
plotted P s (s) for each in Fig. 4. A reasonably good fit to 
all populations is made simply by writing the Gaussian’s 
standard deviation as a = 1/Z. However, the fit is not 
perfect. The structural connectedness of a grain implies 
that a Z-dependency should exist in the cosine factor. 



FIG. 3: Distribution of load anisotropy s depicting single 
grain states. Solid curve — DEM. Dashed curve — theoretical 
prediction, the numerical solution of the transport equation. 
The excellent fit helps validate the theory. 



FIG. 4: The same load anisotropy distribution from the DEM 
as in Fig. 3 except segregated into three graphs by grain coor- 
dination number Z (shifted vertically for clarity). An empiri- 
cal fit was suggested by the theory (dashed curves) which fits 
all three Z populations when the standard deviation cr = 1/Z 
is the only parameter. Despite the disorder of a granular 
packing, this result hints at an underlying elegance in the 
DOS. 


We found from the DEM data that the distribution 
of normalized total loads, Pt(t), is very sensitive to Z 
and does not compare well with the theory in the re- 
gion t < 1. Segregating by Z populations as shown in 
Fig. 5, we see that the error in the region t < 1 is pri- 
marily due to the Z = 3 population, the other popula- 
tions being reasonably well predicted the theory in that 
region. Knowing that the theory’s numerical solution is 
fitted by the functional form P t (t) = t a e~^, we tried 
the same form for all three Z populations, too. We find 
that all three are described excellently by (d = 2Z — 4 
and a = (3 — 1. This can be written neatly by defining 
a convolution operator C m such that C m [/(£)] convolves 
m functions f(t) together. Then P t (t) = C m [e“*] where 
m = 2Z — 4 for all Z. This unifying expression along 
with the beauty of Fig. 5 speaks that this is not an em- 
pirical coincidence. It strongly implies the existence of a 
basis change in the phase space such that the packing’s 
contact forces are resolved into the loading modes of the 
grains, modes which are statistically independent at the 
grain (hence the convolutions) and canonical (hence the 
Gibbs distributions e“*). The number of a grain’s load- 
ing modes m is not generally the number of its contacts 
Z, so the information contained in the new bases is not 
localized in the fabric. However, (m) = (Z) exactly, so 
the change in basis preserves the overall dimensionality 
of the space. This forms an interesting analogy to the 
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FIG. 5: Distribution of total load t depicting single grain 
states in the DEM, segregated into three graphs by their co- 
ordination number Z (shifted vertically for clarity). The the- 
ory does not address the differences between Z populations, 
but it does imply a particular functional form. The empiri- 
cal fits based on that form are Pt(t) = f( 2Z - 5 ) e -( 2Z - 4 ) t > As 
discussed in the text, this implies an elegant result that has 
far-reaching implications. 

molecular vibrations in a solid, which are resolvable into 
non-localized phonons. While we are accustomed to find- 
ing such elegant results in ergodic thermal systems, we 
did not expect to find the elegance of canonical statistics 
lying hidden within the untidiness of a granular packing. 

In summary, it is clear that the transport equation 
needs to be generalized beyond the constant Z = 4 ap- 
proximation, but nevertheless its quantitative predictions 
for Pf(f ), P x (/ x ), and P s (s) are quite good. The predic- 
tion of P t (t) is good, too, when restricted to the relevant 
Z — 4 case. By finding functional fits to the theory we 
discover simple generalizations that describe each of the 
Z populations distinctly, demonstrating simplicity in the 
statistics. An unintended benefit is that this pointed to- 
ward a previously unrecognized canonical picture in the t 
distributions. These are all outstanding confirmations of 
the theory. Returning to the original goal of this Letter, 
we conclude that the Pf(f) seen in a dynamic simula- 
tion is not indicating any dynamical bias in the measure, 
but only the features that Edwards’ hypothesis predicts. 
Therefore, the important granular phenomena that are 
known to produce those features, such as the formation 
of force chains and the diminishing number of very weak 
forces after the jamming transition, are also contained as 
inherent elements of Edwards’ ensemble without recourse 
to the individual grains’ dynamics. Since it is the state 
of these contact forces that determines which grains tip 
and roll during quasi-static rheology, this implies that 
Edwards’ hypothesis may eventually prove sufficient for 
the development of a first-principles rheological theory. 
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